Investigating the modification of the DRHODZ variable in the MSIS.f90 subroutine.¶
DRHODZ is the part of the calculation of the partial derivative of mass density in the D matrix of the variational equations.
This notebook looks at how updating the calculation impacts the DRHODZ term.
Want to show:¶
MSIS00
original DRHODZ
updated DRHODZ
MSIS2
original DRHODZ
updated DRHODZ
Import base parameters¶
[1]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
from plotly.subplots import make_subplots
import plotly.express as px
import copy
import sys
import numpy as np
[2]:
%load_ext autoreload
%autoreload 2
import sys
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/')
from pygeodyn_ReadOutput import pyGeodyn_Readers
run_params = {}
run_params['run_ID'] = 'RUN'
run_params['satellite'] = 'starlette'
run_params['den_model'] = None
run_params['empirical_accels'] = False
run_params['SpecialRun_name'] = None
run_params['options_in'] = {'DRHODZ_update':True}
run_params['verbose'] = False
run_params['arc'] = [ '030914_2wk',
'030928_2wk',
'031012_2wk',
'031026_2wk',
'031109_2wk',
'031123_2wk',
'031207_2wk'
]
Load different versions¶
[3]:
#---MSIS2-DRHODZ-ON---------------------------------
params_msis2_DrhodzOn = copy.deepcopy(run_params)
params_msis2_DrhodzOn['den_model'] = 'msis2'
params_msis2_DrhodzOn['SpecialRun_name'] = ''
params_msis2_DrhodzOn['options_in'] = {'DRHODZ_update':True}
Obj_msis2_DrhodzOn = pyGeodyn_Readers( params_msis2_DrhodzOn )
# Obj_msis2_DrhodzOn.getData_OnAllArcs()
#---MSIS2-DRHODZ-OFF---------------------------------
params_msis2_DrhodzOff = copy.deepcopy(run_params)
params_msis2_DrhodzOff['den_model'] = 'msis2'
params_msis2_DrhodzOff['SpecialRun_name'] = '_drhodzOrig'
params_msis2_DrhodzOff['options_in'] = {'DRHODZ_update':False}
Obj_msis2_DrhodzOff = pyGeodyn_Readers( params_msis2_DrhodzOff )
# Obj_msis2_DrhodzOff.getData_OnAllArcs()
#---MSIS00-DRHODZ-ON---------------------------------
params_msis00_DrhodzOn = copy.deepcopy(run_params)
params_msis00_DrhodzOn['den_model'] = 'msis00'
params_msis00_DrhodzOn['SpecialRun_name'] = ''
params_msis00_DrhodzOn['options_in'] = {'DRHODZ_update':True}
Obj_msis00_DrhodzOn = pyGeodyn_Readers( params_msis00_DrhodzOn )
# Obj_msis00_DrhodzOn.getData_OnAllArcs()
#---MSIS00-DRHODZ-OFF---------------------------------
params_msis00_DrhodzOff = copy.deepcopy(run_params)
params_msis00_DrhodzOff['den_model'] = 'msis00'
params_msis00_DrhodzOff['SpecialRun_name'] = '_drhodzOrig'
params_msis00_DrhodzOff['options_in'] = {'DRHODZ_update':False}
Obj_msis00_DrhodzOff = pyGeodyn_Readers( params_msis00_DrhodzOff )
# Obj_msis00_DrhodzOff.getData_OnAllArcs()
Calling pygeodyn with multiple arcs...
Calling pygeodyn with multiple arcs...
Calling pygeodyn with multiple arcs...
Calling pygeodyn with multiple arcs...
Read Data into Objects¶
[4]:
Obj_msis2_DrhodzOn.getData_OnAllArcs()
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st030914_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st030928_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031012_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031026_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031109_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031123_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff/.../st031207_2wk.goco05s
[5]:
Obj_msis2_DrhodzOff.getData_OnAllArcs()
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st030914_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st030928_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st031012_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st031026_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st031109_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st031123_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis2/msis2_acceloff_drhodzOrig/.../st031207_2wk.goco05s
[6]:
Obj_msis00_DrhodzOn.getData_OnAllArcs()
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st030914_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st030928_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031012_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031026_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031109_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031123_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff/.../st031207_2wk.goco05s
[7]:
Obj_msis00_DrhodzOff.getData_OnAllArcs()
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st030914_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st030928_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st031012_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st031026_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st031109_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st031123_2wk.goco05s
File path:
Loading /data/data_geodyn/results/st/msis00/msis00_acceloff_drhodzOrig/.../st031207_2wk.goco05s
Plot Results¶
[8]:
# Obj_msis00_DrhodzOff.Density['']
# Obj_msis00_DrhodzOff.__dict__['Density']
# X = obj_m1.__dict__['Density'][arc]['Date'][::decimate]
# Y = obj_m1.__dict__['Density'][arc]['drhodz (kg/m**3/m)'][::decimate]
[9]:
%load_ext autoreload
%autoreload 2
import sys
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/')
from pygeodyn_Analysis import plot_drhodz
fig = make_subplots(rows=2, cols=1 )
fig = plot_drhodz(fig, Obj_msis00_DrhodzOff,Obj_msis00_DrhodzOn, 0, 200, 'old drhodz','new drhodz' )
# fig = plot_drhodz(fig, , 1, 200, 'DRHODZ update On')
iplot(fig)
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
[10]:
%load_ext autoreload
%autoreload 2
import sys
sys.path.insert(0, '/data/geodyn_proj/pygeodyn/')
from pygeodyn_Analysis import plot_drhodz
fig = make_subplots(rows=2, cols=1 )
fig = plot_drhodz(fig, Obj_msis2_DrhodzOff,Obj_msis2_DrhodzOn, 0, 200, 'old drhodz','new drhodz' )
iplot(fig)
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
[ ]:
[ ]:
[ ]:
[11]:
# # params_msis00_origDrhoDz
# params_msis2_origDrhoDz
[12]:
# from pygeodyn_ReaderDevelop import pyGeodyn_Readers
# Obj_msis00_origDrhoDz = pyGeodyn_Readers( params_msis00_origDrhoDz )
# Data_msis00_origDrhoDz = Obj_msis00_origDrhoDz.getData_UserChoose(['Density'])
# # Obj_msis2_origDrhoDz = pyGeodyn_Readers( params_msis2_origDrhoDz )
# # Data_msis2_origDrhoDz = Obj_msis2_origDrhoDz.getData_UserChoose(['Density'])
# Obj_msis00_updateDrhoDz = pyGeodyn_Readers( params_msis00_updateDrhoDz )
# Data_msis00_updateDrhoDz = Obj_msis00_updateDrhoDz.getData_UserChoose(['Density'])
# # Obj_msis2_updateDrhoDz = pyGeodyn_Readers( params_msis2_updateDrhoDz )
# # Data_msis2_updateDrhoDz = Obj_msis2_updateDrhoDz.getData_UserChoose(['Density'])
[13]:
# import plotly.graph_objects as go
# import numpy as np
# from plotly.offline import plot, iplot
# # %matplotlib inline
# from plotly.subplots import make_subplots
# import plotly.express as px
# # import plotly
# # plotly.offline.init_notebook_mode(connected=True)
# col1 = px.colors.qualitative.Plotly[0]
# col2 = px.colors.qualitative.Plotly[1]
# col3 = px.colors.qualitative.Plotly[2]
# col4 = px.colors.qualitative.Plotly[3]
# data_nums_1 = 1000
# data_nums_2 = 100
# fig = make_subplots(rows=2, cols=1,
# subplot_titles=("Direct Comparison","Percent Difference in Original and Edited DRHODZ"),
# )
# Obj00 = Data_msis00_origDrhoDz.Density
# Obj01 = Data_msis00_updateDrhoDz.Density
# # Obj10 = Data_msis2_origDrhoDz.Density
# # Obj11 = Data_msis2_updateDrhoDz.Density
# # for arcs in arc_list:
# fig.add_trace(go.Scatter(x=Obj00['Date'][::data_nums_2],
# y=np.abs(Obj00['drhodz (kg/m**3/m)'][::data_nums_2]),
# # name=name_m1,
# mode='markers',
# marker=dict(color=col1,
# size=4,
# ),
# ),
# row=1, col=1,
# )
# fig.add_trace(go.Scatter(x=Obj01['Date'][::data_nums_2],
# y=np.abs(Obj01['drhodz (kg/m**3/m)'][::data_nums_2]),
# # name=name_m1,
# mode='markers',
# marker=dict(color=col2,
# size=4,
# ),
# ),
# row=1, col=1,
# )
# # fig.add_trace(go.Scatter(x=Obj10['Date'][::data_nums_2],
# # y=np.abs(Obj10['drhodz (kg/m**3/m)'][::data_nums_2]),
# # # name=name_m1,
# # mode='markers',
# # marker=dict(color=col3,
# # size=4,
# # ),
# # ),
# # row=1, col=1,
# # )
# # fig.add_trace(go.Scatter(x=Obj11['Date'][::data_nums_2],
# # y=np.abs(Obj11['drhodz (kg/m**3/m)'][::data_nums_2]),
# # # name=name_m1,
# # mode='markers',
# # marker=dict(color=col4,
# # size=4,
# # ),
# # ),
# # row=1, col=1,
# # )
# A0 = (Obj00['drhodz (kg/m**3/m)'][::data_nums_2])
# B0 = (Obj01['drhodz (kg/m**3/m)'][::data_nums_2])
# C0 = ((B0-A0)/A0)*100
# fig.add_trace(go.Scatter(x=Obj00['Date'][::data_nums_2],
# y=C0,
# name = '% diff',
# mode='markers',
# marker=dict(color='green',
# size=5,
# ),
# # showlegend=False,
# ),
# row=2, col=1,
# )
# # A1 = (Obj10['drhodz (kg/m**3/m)'][::data_nums_2])
# # B1 = (Obj11['drhodz (kg/m**3/m)'][::data_nums_2])
# # C1 = ((B-A)/A)*100
# # fig.add_trace(go.Scatter(x=Obj10['Date'][::data_nums_2],
# # y=C,
# # name = '% diff',
# # mode='markers',
# # marker=dict(color='red',
# # size=3,
# # ),
# # # showlegend=False,
# # ),
# # row=2, col=1,
# # )
# fig.update_layout(
# title="Modifying DRHODZ-- Comparison",
# )
# fig.update_yaxes(type="log", exponentformat= 'power', row=1, col=1)
# # fig.update_yaxes( exponentformat= 'power', row=2 , col=1)
# fig.update_yaxes( exponentformat= 'power',)
# fig.update_layout(legend= {'itemsizing': 'constant'})
# fig.update_yaxes(title_text="kg/m^4", row=1, col=1)
# fig.update_yaxes(title_text="%", row=2, col=1)
# fig.update_xaxes(title_text="Date", row=1, col=1)
# fig.update_xaxes(title_text="Date", row=2, col=1)
# fig.update_layout(legend= {'itemsizing': 'constant'})
# # fig.update_xaxes(tickangle=45)
# fig.update_layout(
# autosize=False,
# width=900,
# height=700,
# )
# fig.update_layout(
# font=dict(
# # family="Courier New, monospace",
# size=14,
# ),)
# iplot(fig)
[14]:
# import os
# # if not os.path.exists("msis_update_images"):
# # os.mkdir("msis_update_images")
# fig.write_image("plot_drhodz_msis00_eval.png")
[ ]:
[ ]:
[ ]:
[ ]: